Math. Model. Nat. Phenom. 
2012 



o : Non-homogeneous random walks and subdiffusive transport of cells 

> 
o 

in 



o 

(N 



S. Fedotovlll, A. O. Ivanov " and A. Y. Zubarev " 

" School of Mathematics, The University of Manchester, Manchester, M13 9PL, UK 
Department of Mathematical Physics, Ural Federal University, Ekaterinburg, 620083, Russia 



g . Abstract. This paper is concerned with a non-homogeneous in space and non-local in time 

random walk model for anomalous subdiffusive transport of cells. Starting with a Markov model 

^ I involving a structured probability density function, we derive the non-local in time master equation 

^ ■ and fractional equation for the probability of cell position. We show the structural instability of 

■ fractional subdiffusive equation with respect to the partial variations of anomalous exponent. We 
^ . find the criteria under which the anomalous aggregation of cells takes place in the semi-infinite 

^ ! domain. 
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^ : 1. Introduction 

(N 

■ Random cell movement plays a very important role in embryonic morphogenesis, wound heal- 
ing, cancer cells proliferation, and many other physiological and pathological processes [[34|. The 
microscopic theory of the migration of cells and bacteria towards a favorable environment (chemo- 
taxis) is based on random walk models |[8lll9l[3Tl[32l|. The "velocity -jump" models concern with 

, self-propelled motion involving the runs and tumbles, while "space-jump"models deal with the 

I cells making jumps in space. Much of the literature on the theoretical studies of cells motility 

^ ■ has been concerned with Markov random walk models (see, for example, flU [8l). However, the 

analysis of random movement of wild-type and mutated epithelial cells shows the anomalous dy- 
namics of cell migration [7] (see also IjlTlD . Over the past few years there have been several 
attempts to model non-Markovian anomalous cell transport involving subdiffusion and superdiffu- 
sion |I71|9l[T0l[l2l[T3l[T7]|- In this paper we shall deal with a non-Markovian "space-jump"model 
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that describes the non-homogeneous in space subdiffusive transport of cells. 

First let us consider a Markov model for random cell movement along one-dimensional lattice 
such that all steps are of equal length 1. We define the probability 

p{k,t) = FT{X{t) = k} (1.1) 

that the position of cell X(t) is at point k E Z at time t.We introduce at each point k the rate 
of jump to the left and the rate of jump to the right X{k). This random walk is called a 
generalized birth-death process (H. The master equation for p{k, t) can be written as 

= x{k - l)p{k - 1, t) + fi{k + l)p{k + 1, t) - {X{k) + I2{k))p{k, t). (1.2) 

crt 

This model corresponds to the case when intervals between jumps at point k are exponentially 
distributed with parameter \{k)+n{k) . When the cell makes a jump from the position k, it jumps to 
the right with probability X{k) / {X{k)+ fi{k)) and to the left with the probability fi{k) / {X{k)+ fi{k)) 

m. 

The dependence of /i(A;) and A(A;) on space can be introduced in different ways depending on 
how cells sense the surrounding environment. For the local chemotaxis models, the rates X{S{k)) 
and fj.{S{k)) are the functions of the local concentration of the chemotactic substance S{k). In 
the continuous limit, the master equation (11.21) can be reduced to the classical advection-diffusion 
equation in which the cell flux involves the standard diffusion term and the advection term due to 
chemotaxis. Note that there exist several non-local and barrier models that are different in terms 
of the dependence of rate functions on the chemotactic substance. For example, the rates fi{k) and 
X{k) can depend on the concentration of the chemotactic substance at neighbouring positions k — 1 
and k + 1 (see the details in [|32l l4l). If we consider only positive values of k, we need to implement 
boundary conditions at the point k = 1. Here we assume that if cell hits the wall on the boundary, 
it is reflected with the probability 1 — x and absorbed by wall with the probability x- Then one can 
write At) = (l-A(l)At-/i(l)At)p(l,t)+ /i(l)(l-x)p(l,t)At+/i(2)p(2,t)At+o(At). 

In the limit At — we obtain 

= -xMlMl,^) +M2M2,t) - A(l)p(l,t), (1.3) 

where < x < 1- 

Non-uniform stationary solution of master equation (11.21) can be interpreted as cell aggregation 
phenomenon ||32]| . In particular, if there is no absorption on the boundary (x = 0), the stationary 
solution pst{k) can be easily found from (11.21 ) and (11.31 ). We obtain 

Pstik) = Pst{l)li .^/P^. , k>l, (1.4) 



where 



oo fe— 1 



-1 



k=2 
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provided the series is convergent. 

It is tempting to generalize the master equation (11.21 ) for the anomalous case by replacing the 
time derivative with the Caputo derivative ^ [24]| 

d''p{k,t) ^ 1 f dp{k,u) du 

df" r(i - io du {t - uy-'^ 

as it is done in ll33l for a fractional linear birth-death process. Here u is the anomalous exponent: 
< < 1. Although this generalization is very attractive from a mathematical point of view, it 
is not appropriate for a non-homogeneous medium for which the exponent u depends on k. The 
non-homogeneous fractional equation for p{k, t) can be written as 

^-^^ = a{k-l)Vl~''^'-'^p{k-l,t)+h{k+^^^^^^ 

(1.6) 

where V] "^^"^ is the Riemann-Liouville fractional derivative with varying order 



Here u{k) is the anomalous exponent corresponding to the site k and the anomalous rate coef- 
ficients a{k) and h{k) have to be determined (see (13.101) ). The crucial point here is that the 
anomalous exponent u{k) depends on the site k. The fractional equation (11.61) cannot be rewritten 
in terms of Caputo derivative (11.51) . It turns out that even small non-homogeneous variations of 
the exponent u lead to a drastic change of p{k,t) in the limit t — oo lfT4l . It means that the 
subdiffusive fractional equations with constant anomalous exponent u are not structurally stable. 
If, for example, the point k = M has the property that z/(Af) < i'(k) for all A; 7^ M and x = 0, 
one can find that 

p{k, t) -> 0, p{M, t) ^ 1, l<k<N (1.8) 

as t — 00. This result has been interpreted as anomalous aggregation of cells at the point k = M 
lfT2l . In this paper we shall find the conditions for anomalous aggregation for the semi-infinite 
interval 1 < /c < 00. It should also be noted that non-homogeneous variations of the exponent 
v destroy the Gibbs-Boltzmann distribution as a long time limit of the fractional Fokker-Planck 
equation |T4|. 

Another extension of traditional Markov random walks models is non-Markovian theory of 
anomalous transport with reaction dynamics lITTl l28l [30ll35l l36l |20| . In particular, this theory has 
been used for the analysis of the proliferation and migration dichotomy of cancer cells [|9l [T0l[T3l . 
In this paper we consider the inhibition of cell growth by anticancer therapeutic agents. To model 
this inhibition we introduce the random death process with non-uniform death rate parameter. We 
assume that during time interval (t, t + At) at point k each cell has a chance 9{k)At + o{At) of 
dying, where 9{k) is the death rate. It is easy to take into account this process for Markov models. 
We just add the term —6{k)p{k,t) to the right hand side of the master equation (11.21) . On the 
contrary, the anomalous master equation involves a non-trivial combination of transport and death 
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kinetic terms because of memory effects IfTSl |28l [Til. In this paper we shall derive the following 
fractional equation 

= a{k - iy-o[k-i)tjyi~^ik^i) ^^(^ _ 1^ tyik-i)t^^ 
+b{k + i)e-e(fe+i)*pi-Kfc+i) + i,t)e^('=+i)*] 
- {a{k) + b{k)) e-eW*r'i-'^('=) [p{k, t)e^(^)*] - e{k)p{k, t). (1.9) 

Instead of the probability p{k,t) for an individual cell one can consider the mean density of 
cells p{x, t) as a function of space x and time t. Then the master equation (11.21) for the density 
p{x, t) can be written as 

^P^^lll = x{x - l)p{x - /, t) + + l)p{x + l,t)- (A(x) + /i(x))p(x, t) - 9{x)p{x, t), (1.10) 

where / is the jump size, 9{x) is the death rate. The advantage of this equation is that one can 
easily take into account various non-linear effects by assuming the dependence of the rate functions 
X{x), p{x) and 6{x) on the average density p{x, t). 

In the anomalous subdiffusive case, the master equation for mean field p{x, t) becomes non- 
local in time. It can be written as a mass balance equation 

= -/(x, t) + I{x -l,t)- 9{x)p{x, t), (1.11) 

where I{x, t) is the total flow of cells from the point xtox + l 

I{x,t) = a{x)e-'^'^^'V]-'^''^ [e'^^^'p{x,t)] - 6(x + /)e-^(^+')*I)^"("+') [e'^^+'^'p{x + l,t)] 

(1.12) 

and I{x — l,t) is the total flow of cells from the point x — I to x 

I{x -l,t) = a{x - /)e-^(--0<i?i-K-0 [e^(--Otp(a; _ i^t)] - 6(x)e-^(^)*I),'-"(") [e"^"^^' p{x , t)] . 

(1.13) 

Here a(x) and b(x) are the anomalous rate functions (see ( 13.101 ) ). One can see that the flow of cells 
I{x, t) depends on the death rate 0{x). It means that in the anomalous case one cannot separate the 
transport of cells from the death process [18|. This phenomenon does not exist in the Markovian 
case. For the Markov model (11.101) the flux /(x, t) is independent from 0{x) : 

/(x, t) = X{x)p{x, t) — + l)p{x + l,t). 

When the density p{x, t) is conserved (6 = 0), the master equation (|1.11|) can be approximated by 
the fractional Fokker- Planck equation |l2l |25l |26l 



dp{x,t) d 
dt dx 



l{a{x)-b{x))Vl-"'''>p{x,t) 



dx"^ 



^\{x) + b{x))Vl-'^''^p{x,t) 



(1.14) 



This is an example of the fractional equation with varying anomalous exponent |l5l. Note that 
a{x) — b(x) ~ / as / — (see (13.161) and (13.181) ). The purpose of the next section is to set up 
a non-Markovian discrete-space random walk model describing cell motility involving memory 
effects, the death process and subdiffusive transport. 
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2. Non-Markovian discrete-space random walk model 

2.1. Random cell motility 

There exist numerous mechanisms that facilitate random cell movement ||34l . In this paper we 
adopt the following random model of cell motility. When the cell makes a jump to position k, the 
time the cell spends here before it makes a jump to point A; — 1 or + 1 is random. It is called the 
residence time or waiting (holding) time. We define the residence time at position k as 

Tk = mm{Ti:,T^), (2.1) 

where and are the independent random times of jump to the left and right respectively. The 
idea here is that there exist internal cellular signals involving two "hidden" independent random 
alarm clocks. If one of the clocks goes off first, say < T^, the cell moves to the right to the point 
k + 1. The other clock "tells" the cell to move left to the point /c - 1 if it goes off first (T^ < T^). 
Note that migration of cells is a highly complicated dynamic process which is regulated by both 
intercellular signals and the surrounding environment. Since we do not know the exact mechanism 
of cell motility we use a stochastic approach involving two random times and for jumping 
to the left and right. Note that if the random times and are exponentially distributed with 
the rates ii{k) and \{k) respectively, we have a classical Markov model with the master equation 
( 11.21) . If the random variables and are not exponentially distributed, the standard Markov 
approach does not work. In this section we consider the non-Markovian case when and are 
independent positive random variables with general survival functions 

^^{k, t) = Pr {Ti: > r} , ^x{k, r) = Pr {T,^ > r} . (2.2) 

The Markov model (11.21) corresponds to the following choice 

^^{k, t) = e-^(^)^ ^xik, t) = e-^^'^)^ (2.3) 

It is convenient to introduce the rate of escape (hazard function) 7(A;, r) from the point k as 

^{k, r) = lim PHr<T,<r + hW,,^} _ 

h^O h 

If we denote the survival function at the point k as 

^(A;,r) = Pr{Tfc >r} 
and the residence time probability density function as 

thenO 
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Now we determine this rate function in terms of statistical characteristics of random residence 
times and T^. It follows from the definition of the residence time at position k (12.11 ) that the 
survival function \I/(/c, r) can be written as a product 

^(/c,r) = ^A(A;,r)^^(A;,r), 

where '^\{k, r) and r) are defined by (12.21) . Differentiation of this equation with respect to 

T gives 

ij{k,T) = ^x{k,T) + ^^{k,T), (2.6) 
where the transition probability density functions il)\{k, r) and il)^{k, r) are defined as 



= (2.7) 



The formula (12.61 ) is the particular case of the general expression for the residence time PDF in 
terms of the transition probability densities (see formula (5) in the classical paper by van Kampen 
tmX ). These transition PDF's have a clear probabilistic meaning. For example, ^/'^(/c, r)Ar is the 
probability that the cell's jump to the left occurs in the time interval (r, r + Ar) since the cell 
arrived at point k. If we divide both sides of (12.61) by the survival function \E'(A;, r) and use the 



formula (12.51) . we obtain 

7(fc,r) = A(A;,r) + /i(/c,r), (2.8) 
where the rate of jump to the right A(A;, r) and the rate of jump to the left r) are defined as 

Xik r) - ^"^^'^^ uik r) - ^^^^-^^ (2 9) 

Note that the transition rates \{k, r) and ii{k, r) can be introduced from the very beginning as it is 
done in [|22ll . By using (12.51) and (12.81) . we write the survival function ^(A;, r) as 

vl>(A;,r) = e"-^o'^{'='^)'^^-/(rMfc.^)'i^_ (2.10) 
The residence time probability density function V^(A;, r) takes the form 



ij{k,T) = (A(A;,r) +/i(A;,r))e-^o^^(^'")'^-/oV(fc,r)dr_ 

For the Markov case for which \{k) and uik) are independent of the residence time variable r, we 
obtain from (12.101) the standard exponential survival function 

corresponding to the Markov master equation (11.21) . 
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2.2. Structured probability density function 

If the residence time probability density function ^ is not exponential, the random walk is non- 
Markovian. The standard method to deal with non-Markovian stochastic processes is to add aux- 
iliary variables to the definition of the random walk to make the process Markovian Here we 
introduce the structured probability density function ^(k, t, r) involving residence time r as auxil- 
iary variable. The structural density gives the probability that the cell position X{t) at time t is at 
the point k and its residence time Tk at point k is in the interval (r, r + dr). This is a standard way 
to deal with non-Markovian random walks [|6l|28]|. Suppose that cells die at random at rate 6{k) 
that depends on k. The density ^{k,t,T) obeys the balance equation 

§ + 1^ = -Kk, r)e - r)e - e{k)C (2.1 1) 

We consider only the case when the residence time of random walker at t = is equal to zero, so 
the initial condition is 

ak,0,T)=poik)6iT), (2.12) 

where poik) = Pr {X{0) = k}. The boundary condition in terms of residence time variable (r = 
0) can be written as ^ 

^{k,t,0)= [ X{k-l,T)^{k-l,t,T)dT+ [ + l,r)^(fc + l,t,r)rfr. (2.13) 
Jo Jo 

In what follows we consider only positive values of k. In which case, we have to specify the 
boundary condition for A; = 1 . We write 

e(l,t,0) = (l-x) / /i(l,r)^(l,t,r)dr+ /" /i(2, r)e(2, t, r)dr. (2.14) 
Jo Jo 

This equation tells us that when cells escape from the point k = 1 and move to the left with the rate 
r), they are adsorbed by the wall with probability x, and reflected back to the position k = 1 
with the probability 1 — x- Note that this boundary condition can be written in many different ways, 
for example, the cells can be reflected to state k = 2. One can also introduce a residence time PDF 
for a wall such that the reflection is not instantaneous. 
We solve (12.111) by the method of characteristics 

^(A;,t,r) = e(A;,t-r,0)e-^o"^{^^'")'^-/oV{fc,r)dr-e(fc)r^ ^ ^ ^ > ^_ ^-2.15) 

The structural density ^ can be rewritten in terms of the survival function r) (|2.10l) and the 
integral arrival rate 

j{k,t)=^ik,t,0) 

as 

^{k,t,T) = j{k,t-T)^{k,T)e-^^''^\ T<t, k>l. (2.16) 
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Our purpose now is to derive the master equation for the probability p{k, t) = Pr {X{t) = k} : 

p{k,t)= [ ^{k,t,T)dT, k>l. (2.17) 
Jo 

Let us introduce the integral escape rate to the right ix{k, t) and the integral escape rate to the left 

if,{k,t) as 

tx{k,t)= [ X{k,T)^{k,t,r)dr, i^{k,t)= [ fi{k,T)^{k,t,T)dT. (2.18) 
Jo Jo 

Then the boundary conditions (12.131) and (12.141) can be written in a very simple form: 

j{k,t) =tx{k-l,t)+tf,{k + l,t), k>l (2.19) 

and 

j(l,t) = (l-x)v(l>^)+V(2,^)- (2.20) 
It follows from ^3, (I02l) . (ITTel) and (I2l8]) that 

^x{k, t) = f T)j{k, t - T)e-'^'^^dT + Mk, t)po{k)e-'^^^\ (2.21) 

i,{k,t) = f^,{k,T)jik,t-T)e-"'''^^dT + Mk,t)po{k)e-'^''^\ (2.22) 
Substitution of (12.121) and (12.161) to (I2.17L gives 

p{k,t)= [ ^{k,T)j{k,t-T)e-^^'''>^dT + ^{k,t)po{k)e-^''''^\ (2.23) 
Jo 

It is convenient to introduce the integral escape rate i{k, t) as the sum of the escape rate to the right 
ix{k, t) and the escape rate to the left i^{k, t) 

i{k,t)=ix{k,t)+i^{k,t). (2.24) 

The balance equation for t) can be written as 

dp{k, t) 



dt 



-i{k,t)+j{k,t)-e{k)p{k,t), k>l. (2.25) 



To obtain a closed equation for p{k, t) we need to express i{k, t) and j{k, t) in terms of p{k, t). 
By applying the Laplace transform il){k, s) = ijj{k, T)e~^'^dT to (12.211) . (I2.22|) and (|2.23l) . we 
obtain 

h{k, s) = i)x{k, s + e{k)) [j{k, s) + po{k)\ , 
i^,{k, s) = ^^{k, s + e{k)) [j{k, s) + po{k)] 
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and 

p{k, s) = i!{k, s + 9{k)) [j{k, s) + po{k)] . 
In the Laplace space we have the following expressions for escape rates 

-^{k,s + e{k)) -^{k,s + e{k)) 

Using inverse Laplace transform and shift theorem we find 

^x{k,t) = f Kxik,t-r)e-''^''^^'-^^p{k,T)dT, 
Jo 

i^.{k,t) = I K,{k,t-T)e~'^'''>^'-^^p{k,T)dT, (2.27) 

where K\{k, t) and K^(k, t) are the memory kernels defined by Laplace transforms 

/f.(fc,.) = «M. k,(k,s)^MM. (2.28) 

It follows from (12.191) . (I2.24L (12.251) and (12.271) that the master equation for the probability p{k, t) 
is 

MM = [' K^{k-l,t-T)p{k-l,T)e-"'''-'^^'-^UT 
ot Jo 

+ [ K^{k + l,t-T)p{k + l,r)e-'^''^'^^'-^^dT 
Jo 

- [ [K,{k, t-T) + K^{k, t - rMk, r)e-^(^')(*-^)(ir - 9{k)p{k, t) (2.29) 

for k > 1. The balance equation for A; = 1 is 

dp{l,t) 



dt 



t) - t) + 1^{2, t) - e{i)p{i, t) 



or 

dt 



Jo Jo 

+ f K,{2, t - r)e-^(2)(i-.)^(2, r)cir - ^(l)p(l, t), (2.30) 
Jo 

where < x < 1- The master equation for p{k, t) can be rewritten in terms of the probability flux 
I{k, t) from the point A; to A; + 1 

I{k,t)= [ Kx{k,t-T)p{k,T)e-^^''^'^'~^UT- [ i^^(A: + l,t-r)p(A; + l,r)e-^('=+^)(*-")c/r 
Jo Jo 

(2.31) 
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as 

= _/(A;, t) + lik - 1, t) - 9ik)pik, t). (2.32) 
In the next section we shall derive fractional master equation for p(k,t). 

3. Anomalous subdiffusion in heterogeneous media 

We now turn to the anomalous subdiffusive case. We assume that the longer cell survives at point 
k, the smaller the transition probability from k becomes. It means that the transition rates \{k, r) 
and n{k, r) are decreasing functions of residence time r. We assume that 

ro(fc) + r To{k) + T 

where To{k) is a parameter with units of time. Both i'x{k) and i'^{k) play a very important role in 
what follows. From (I2.10|) and (|3.1I) we find that the survival function has a power-law dependence 

where the exponent 

u{k) = vx{k) + v^{k) (3.2) 

depends on the state k. Residence time probability density function ^(k,r) = — 9\l/(/c, T)/dT has 
the Pareto form 

The anomalous subdiffusive case |l2l|28]| corresponds to 

iy{k) = iyx{k) + iy^{k) < 1. 

We can notice from (13.11 ) that the ratios X{k, r) and fi{k, r) to X{k, r) + fi{k, r) are independent of 
the residence time variable r that is 

A(/c,r) i^x{k) (J^ik,T) v^{k) 



X{k,T) + fi{k,T) ux{k) + u^{ky X{k,r) + iJ,{k,T) i^x{k) + i/f^ik)' 
In this case it is convenient to introduce the probabilities of jumping to the right 

and to the left 

i^x{k) + v^[k) 
10 
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Note that these jump probabilities are completely determined by the anomalous exponents uxi^k) 
and ^^{k). Note that in the standard CTRW theory, these jump probabilities are given indepen- 
dently [I2l|28l. The corresponding transition PDF's ipx{k,T) = \{k,T)'^(k,T) and ijj^{k,T) = 
fi{k, r)\l/(fc, r) can be rewritten in terms of ^/'(A;, r) 

^xik, t) = px{k)ip{k, r), tjj^ik, r) = p^{k)ip{k, r). (3.6) 

The asymptotic approximation for the Laplace transform of the waiting time density ^{k, r) of the 
Pareto form (|3.3I) can be found from the Tauberian theorem [|T5l 

^p{k,s) c^l- g{k)s''^''\ s^O 

with 

g{k) = T{l-u{k)){To{k)r^'K (3.7) 
We obtain from (12.281) and (13.61 ) the Laplace transforms of the memory kernels 

Kx{k,s)^^^^^^^^, .^0. (3.8) 

Therefore, the integral escape rates to the right ix{k, t) and to the left i^(/c, t) in the subdiffusive 
case are 

xx{k^t) = a{k)e-'('^'v'r^'^[pik,t)e'('^'], 

i^{k,t) = 6(A;)e-^«*P^"(') [p(A;,t)e^('=)*] . (3.9) 

Here d]'"^^^ is the Riemann-Liouville fractional derivative with varying order defined by (11.71) . 
The anomalous rate functions a{k) and h{k) are 

a(k) = ,, , = 7TT, (3.10) 

9{k) u{k)Til-uik)){Toik)Y^''^ 

Uk) = pM = (3 11) 

9{k) iy{k)T{l-iy{k)){To{k)y^''^ 

with the anomalous exponent ^{k) defined in (13. 21) . The master equation (12.291) takes the form of 
non-homogeneous fractional equation 



dp(k, t) 

= a 



dt 

+b{k + l)e~'^'+'^'p{k + l)pi-^('=+^) [p{k + l,t)e^('=+i)*] 

- {a{k) + b{k)) e-'^'^^'vl"'^''^ [p{k, t)e'(*^)*] - e{k)p{k, t) (3.12) 



for k > 1. 
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For A; = 1 with 6{k) = x = 0, we obtain 

= 6(2)P,^-^%(2, t) - a{l)V'r'^'^p{l, t). (3.13) 
The fractional probability flux Iy{k, t) from the point A; to A; + 1 is 

(3.14) 

The equation (13.121) can be rewritten in terms of the probability flux I,^{k, t) as 

dpik, t) ^ ^ ^^^^ - 1, t) - e{k)p{k, t). (3.15) 



dt 



Note that in the continuous case (x instead of k), the anomalous advection l(a(x) — b{x)) in 
the fractional Fokker-Planck equation (|1.14l) is proportional to the difference in the anomalous 
exponents i'x{x) and u^i^x), since 

aix) - b(x) = P'^"^^ - P^^""^ = '^^(^^ - '^^'(^^ (3 16) 

9{x) u{x)T{l - u{x)) (ro(x))"(") 

The difference i^xi^x) — i^^{x) can be approximated in the different ways. For example, in the 
case of chemotaxis this difference is proportional to the gradient of the local concentration of 
the chemotactic substance S{x). Let us show this. We consider the standard non-local model 
for which the jump probabilities px{k) and p^{k) depends on the chemotactic substance S{k) as 
follows AM II 

p,{k) = Ae^C^e^-D-^W), (3.17) 
where the parameter A is determined from px{k) + p^{k) = 1. Then [fTTll 



Px{k) -P,,{k) 



g/35(fc+l) _ g/35(fc-l) 
^pS(k+l) g^S(fc-l) • 

In the continuous case 

In the limit / — t- 0, we obtain the standard chemotaxis model 



(99 

Px{x)-p,{x) = pi— + o{l) (3.18) 



for which the chemotaxis equation takes the form 



dp{x,t) d 
dt dx 



g{x) dx 



dx^ 



2 



(3.19) 



In the next section we will be concerned with the asymptotic behavior of p(k, t) as t — )■ 00 for 

e{k) = 0. 
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4. Structural instability and anomalous aggregation 

It has recently been shown that the subdiffusive fractional equations with constant anomalous expo- 
nent u ma bounded domain [0, L] are not structurally stable with respect to the non-homogeneous 
variations of parameter u lfT4]| . It turns out that the spatial variations of the anomalous exponent 
lead to a drastic change in asymptotic behavior of p{k, t) for large t. The purpose of this section is 
to find the conditions of this structural instability in semi-infinite domain 1 < A; < oo. We consider 
the case when 6{k) = x = for which the total probability is conserved 

oo 

Y,P{k,t) = l (4.1) 

k=l 

and the fractional probability flux Iy{k, t) from the point /c to A; + 1 is 

h{k, t) = a{k)Vl~''^^^p{k, t) - b{k + l)p{k + l)Vl^'^''^'^p{k + l,t). (4.2) 

For simplicity we assume that the initial conditions are j9o(l) = 1 and po{k) = for k ^ 1. Taking 
the Laplace transform of (11.61 ) and ( 14.11 ) we obtain 

sp{k,s) = a{k-l)s^-'''^''-^'>p{k-l,s) + b{k + l)s^-''^''+^^p{k + l,s) 

-{a{k) + h{k))s'-''^^^p{k,s) (4.3) 

and 

oo 

Y,sp{Ks) = l. (4.4) 

fc=i 

Since there is no flux of cells outside the left border, we have for /c = 1 

spil, s)-l= 6(2)5^-^(2)^(2, s) - a(l)s^-'^Wp(l, s). (4.5) 
In the limit s — )■ 0, one can obtain from ( 14.51) simple formula expressing p(2, s) in terms of ]5(1, s) 

p(2, s) ~ -^^J^^ s), s ^ 0. 

In general, we find from (14.31) and (14.51 ) p{k, s) in terms of p{k — 1, s) 

p{k,s)c^^ p{k-l,s), k>l, s^O. (4.6) 

b{k) 

This formula has very simple probabilistic meaning: the flux I^{k — 1, t) — )■ as t — )■ oo. If we 
take the Laplace transform of I^{k — l,t) from ( 14.21) . we obtain 

t{k -l,s) = a{k - l)si-'^('=-i)p(A; - 1, s) - h{k)p{k)s^-''^^^p{k, s). (4.7) 

It follows from ( |431 ) that iy{k, s) ~ as s 0. 
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4.1. Stationary solution to fractional equation with constant anomalous ex- 
ponent 

Let us assume that the anomalous exponent z/(A;) is independent of the position k that is z/ = const. 
Let us find stationary solution to the fractional master equation (|1.6I) 



Pst{k) = lira p{k,t) = Urn sp{k, s). (4.8) 



It follows from (|43I) that 



a(k — 1) 

p{k,s) - \ ^ p{k-l,s), k>l, s^O. 
b{k) 



or 

k-l . .X 

^^j^ + 1) 

Using the normalization condition (14.41) and (14.81) . we obtain the stationary solution of the equation 

fJl Kj + 1) 

where 

(OO k—1 / .N 

A;=2 

If the sum 



oo fc— 1 / .X 

is divergent, the stationary solution does not exist. In particular, if we assume that the anomalous 
rate functions a and b are equal, that is, a{k) = b{k + 1), then for a finite domain with k = 
1, 2, N, we obtain uniform distribution p^t (A;) = 1 /N for every k. The stationary solution (14.101) 
is very similar to (11.41) corresponding to the Markov birth-death model. However, this similarity 
is very deceptive, because (14.101) is not structurally stable with respect to the non-homogeneous 
variations of parameter z/. The aim of next subsection is to show this structural instability. 

4.2. Anomalous aggregation. 

Now we consider non-homogeneous case for which the anomalous exponent depends on k. We 
assume that the point k = M has the property that i^(Af) < z/(A;) for all k ^ M. Our purpose now 
is to find the conditions under which 

limp(M,t) = l, limp(fc,t)=0, k M. (4.11) 

t—^oo t—^oo 



14 



S. Fedotov et al. 



Non-homogeneous random walks and anomalous transport of cells 



It means that the total probability concentrates just at one point k = M. This phenomenon is 
called an anomalous aggregation [|T2ll . This asymptotic behavior of cells was observed in exper- 
iments on phagotrophic protists when "cells become immobile in attractive patches, which will 
then eventually trap all cells" llT6l . In the Laplace space, (14.1 II) takes the form 

\imsp(M, s) = 1, \imsp(k, s) = 0, k^M. 
We can rewrite the normalization condition (14.41) as 

M-l oo 

sp{M, s) + ^ sp{M -k,s) + ^ sp{M + k,s) = 1. (4.12) 

k=l k=l 

By using ( 14.61) . we express p(M + k, s) in terms of p{M, s) as follows 

k 

p{M + fc, g) - p{M, s) II "^^J: ^ T gKM+fc)-.(M)^ ^_^0^ (4_j3) 

j=i ^ 

Now we write the formula for p{M — k, s) in terms of p{M, s) 

k 

p{M -k,s)c. piM, H ^^^,~/\^\ k = l,...,M-l, 

(4.14) 

Now we substitute (14.131) and (14.141) into (14.121) and use sp{M, s) as a common factor 

sv(M s)(l + V A KM -3 + 1) ^ u{M+k)~u(M) A a(M + j - 1) \ ^ 

Since i^(A'f) < u{k) for any ky^M,wo have s'^(^^+fc)-KA^) ^ q and s'^(A^-fc)-KA^) ^ q as s ^ 0. 
We conclude that if 

\^ u{M+k)-u{M) TT q(M + j - 1) 

as s — > 0, then sp{M, s) — j- 1, while sp{k, s) — 0. It means that in the limit t — > oo, we obtain 
(14.1 II) . If instead of probability p{k, t) we consider the mean density of cells p (x, t) , the formula 
(14.1 II) can be rewritten as p (x, t) — > 5{x — Xmm) as t -> oo, where Xmin is the point on the interval 
[0, oo) at which i/{x) takes its minimum value. Note that this result was obtained for a symmetrical 
random walk in the context of chemotaxis and anomalous aggregation [[T2|. 

5. Conclusions. 

We have studied a non-homogeneous in space and non-local in time random walk model describ- 
ing anomalous subdiffusive transport of cells. Using a Markov model with structured probability 
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density function, we have derived non-local in time and fractional master equations for the proba- 
bility of cell position. The advantage of our probabilistic approach is that it allows us to take into 
account the death process within the general non-Markovian random walk. The main feature of 
our fractional model is that the transition probabilities for jumping on the left and right depend 
inversely on the residence time variable. This dependence induces power-law residence time dis- 
tribution and ultimately the anomalous subdiffusion of cells. It has recently been shown that the 
subdiffusive fractional equations with constant anomalous exponent are not structurally stable in a 
bounded domain with respect to the non-homogeneous variations. In this paper we have extended 
and complemented our previous results for infinite domain and found exact conditions under which 
the structural instability takes place. Our model can be generalized in many ways, e.g., by mod- 
elling the residence time by internal chemical reactions via a stochastic or ordinary differential 
equations instead of simple equation for the residence time dr/dt = 1. It would be interesting to 
take into account the density-dependent dispersal ||29ll including non-linear exclusion process with 
cell-to-cell adhesion (HI 121. 
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